
#include <msParser.h>
#include <msZMat.h>
#include <msKineticOperator.h>
#include <msScannerPes.h>
#include <msSampler.h>
#include <eigen3/Eigen/src/Geometry/Translation.h>
#include <qcCfour.h>
#include <time.h>

using namespace atomism;
  
class msSamplerN4 : public atomism::msSampler
{ 
public:  
     msSamplerN4() : msSampler() { constructVar("msSamplerN4","SamplerN4","sampler N4"); }  
     
     void update()  { msSampler::update(); 
                      };		     
     void initialize() {  msSampler::initialize();  
                        };			    
     static boost::shared_ptr<msSamplerN4> New(){ boost::shared_ptr<msSamplerN4> T( new msSamplerN4 );
						  T->initialize(); T->update(); return T;
						 }     
     bool getNextPoint(boost::shared_ptr<msGeneralizedCoordinates> coor){
       
              if( Scanner->getLastEnergy()<500  )	      
		(*coor)[0]->Value-= max(0.05 + ((*coor)[0]->Value-2)/10.,0.05);
	      else{ coor->setRandom();(*coor)[0]->Value=6.;
	      }
	      if(nScanned==0){
		(*coor)[0]->Value=6;(*coor)[1]->Value=0.1;
		(*coor)[2]->Value=0.1;(*coor)[3]->Value=0.1;
	      }
	      KineticOperator->setDynamicDofs();
	      ZMat->updatePositions(); 
	      nScanned++;
	      return 1;
              };
  boost::shared_ptr<msZMat> ZMat;
  boost::shared_ptr<msKineticOperator> KineticOperator;  
  boost::shared_ptr<msScannerPes> Scanner;
};

int main(int argc, char **argv) {
 
 MPI_Init(NULL, NULL);
 int world_rank; int world_size; int name_len;
 char processor_name[MPI_MAX_PROCESSOR_NAME];
 MPI_Comm_rank(MPI_COMM_WORLD, &world_rank); 
 MPI_Comm_size(MPI_COMM_WORLD, &world_size);
 MPI_Get_processor_name(processor_name, &name_len);
 
 using namespace atomism;

 boost::shared_ptr<msUnitsManager> unitMng 	= 		msUnitsManager::New();

 unitMng->set( "Angstrom Degree amu hartree ps" );

 boost::shared_ptr<msZMat> zmt			=  		msZMat::New(unitMng);
 zmt->setId("zmat");
 zmt->set("X N 0 0.55 N 0 0.55 1 179.9 X 0 5. 1 90. 2 0. N 3 0.55 0 90. 1 180. N 3 0.55 4 179.9 0 180.");
 
 boost::shared_ptr<msSystem> system =  		msSystem::New(unitMng);
 system->addEntity(zmt);

 boost::shared_ptr<msScalarVariable> q0 = msScalarVariable::New("Angstrom");
 q0->set(5,0.5,6,0.1,0.1); q0->setId("q0");
 boost::shared_ptr<msScalarVariable> q1 = msScalarVariable::New("Degree");
 q1->set(120,0.1,179.9,2,2); q1->setId("q1");
 boost::shared_ptr<msScalarVariable> q2 = msScalarVariable::New("Degree");
 q2->set(120,0.1,179.9,2,2); q2->setId("q2");
 boost::shared_ptr<msScalarVariable> q3 = msScalarVariable::New("Degree");
 q3->set(120,0.,360,2,2); q3->setId("q3");

 
 boost::shared_ptr<msGeneralizedCoordinates> gcoors = 
					     msGeneralizedCoordinates::New(unitMng);
 gcoors->addVar(q0); gcoors->addVar(q1); gcoors->addVar(q2); gcoors->addVar(q3);

 boost::shared_ptr<msParser> kinfct0 = 	msParser::New("Angstrom");
 kinfct0->setCoordinates(gcoors);kinfct0->setExpression("q0");kinfct0->setId("f=q0");
 boost::shared_ptr<msParser> kinfct1 = 	msParser::New("Degree");
 kinfct1->setCoordinates(gcoors);kinfct1->setExpression("q1");kinfct1->setId("f=q1");
 boost::shared_ptr<msParser> kinfct2 = 	msParser::New("Degree");
 kinfct2->setCoordinates(gcoors);kinfct2->setExpression("q2");kinfct2->setId("f=q2");
 boost::shared_ptr<msParser> kinfct3 = 	msParser::New("Degree");
 kinfct3->setCoordinates(gcoors);kinfct3->setExpression("q3");kinfct3->setId("f=q3");
 
 boost::shared_ptr<msKineticOperator> kinop =  	msKineticOperator::New(unitMng);
 kinop->set(gcoors,system);
 kinop->addDynamicDof("zmat",3,kinfct0);
 kinop->addDynamicDof("zmat",4,kinfct1);
 kinop->addDynamicDof("zmat",7,kinfct2);
 kinop->addDynamicDof("zmat",8,kinfct3);
 
 boost::shared_ptr<msSamplerN4> sampler  = 	msSamplerN4::New();
 sampler->ZMat=zmt; sampler->KineticOperator=kinop;
 
 boost::shared_ptr<qcCfour> cfour  = 		qcCfour::New();
 cfour->setLevel( qcCfour::G6_31 , qcCfour::CCSDpT );
 cfour->setGenBasSymLink("/org/centers/pecos/APPLICATIONS/cfour/basis/GENBAS");
 cfour->setWorkingDirectory("Cfour_"+output::getString<size_t>(world_rank));
 
 boost::shared_ptr<msZMat> zmtN2   =  		msZMat::New(unitMng);
 zmtN2->setId("zmat"); 
 zmtN2->set("N N 0 1.1");
                      
 double eN2 = -1;
 try{ eN2 = cfour->singlePoint(*zmtN2);
     }
 catch(...){eN2=-1;}
 
 cout<<"Processor "<<world_rank<<": eN2="<<eN2<<endl;
 srand(time(NULL)+2*world_rank);
 time_t begin;time(&begin);
 boost::shared_ptr<msScannerPes> scanner = 	msScannerPes::New();
 scanner->set(zmt,kinop,sampler,cfour);
 scanner->OutputFile="Database_"+output::getString<int>(world_rank);
 scanner->setOrigin(2*eN2);
 sampler->Scanner=scanner;
 scanner->execute(10000);

 time_t end;time(&end);
 cout<<end-begin<<endl;
 MPI_Finalize();
 return 0;
 }


 
